SUBROUTINE svdcmp_sp(a,w,v)

  USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,outerprod
  USE nr, ONLY : pythag
  IMPLICIT NONE

  REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a
  REAL(SP), DIMENSION(:), INTENT(OUT) :: w
  REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v
  INTEGER(I4B) :: i,its,j,k,l,m,n,nm
  REAL(SP) :: anorm,c,f,g,h,s,scale,x,y,z
  REAL(SP), DIMENSION(size(a,1)) :: tempm
  REAL(SP), DIMENSION(size(a,2)) :: rv1,tempn

  m=size(a,1)
  n=assert_eq(size(a,2),size(v,1),size(v,2),size(w),'svdcmp_sp')
  g=0.0
  scale=0.0
  do i=1,n
     l=i+1
     rv1(i)=scale*g
     g=0.0
     scale=0.0
     if (i <= m) then
        scale=sum(abs(a(i:m,i)))
        if (scale /= 0.0) then
           a(i:m,i)=a(i:m,i)/scale
           s=dot_product(a(i:m,i),a(i:m,i))
           f=a(i,i)
           g=-sign(sqrt(s),f)
           h=f*g-s
           a(i,i)=f-g
           tempn(l:n)=matmul(a(i:m,i),a(i:m,l:n))/h
           a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
           a(i:m,i)=scale*a(i:m,i)
        end if
     end if
     w(i)=scale*g
     g=0.0
     scale=0.0
     if ((i <= m) .and. (i /= n)) then
        scale=sum(abs(a(i,l:n)))
        if (scale /= 0.0) then
           a(i,l:n)=a(i,l:n)/scale
           s=dot_product(a(i,l:n),a(i,l:n))
           f=a(i,l)
           g=-sign(sqrt(s),f)
           h=f*g-s
           a(i,l)=f-g
           rv1(l:n)=a(i,l:n)/h
           tempm(l:m)=matmul(a(l:m,l:n),a(i,l:n))
           a(l:m,l:n)=a(l:m,l:n)+outerprod(tempm(l:m),rv1(l:n))
           a(i,l:n)=scale*a(i,l:n)
        end if
     end if
  end do
  anorm=maxval(abs(w)+abs(rv1))
  do i=n,1,-1
     if (i < n) then
        if (g /= 0.0) then
           v(l:n,i)=(a(i,l:n)/a(i,l))/g
           tempn(l:n)=matmul(a(i,l:n),v(l:n,l:n))
           v(l:n,l:n)=v(l:n,l:n)+outerprod(v(l:n,i),tempn(l:n))
        end if
        v(i,l:n)=0.0
        v(l:n,i)=0.0
     end if
     v(i,i)=1.0
     g=rv1(i)
     l=i
  end do
  do i=min(m,n),1,-1
     l=i+1
     g=w(i)
     a(i,l:n)=0.0
     if (g /= 0.0) then
        g=1.0_sp/g
        tempn(l:n)=(matmul(a(l:m,i),a(l:m,l:n))/a(i,i))*g
        a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
        a(i:m,i)=a(i:m,i)*g
     else
        a(i:m,i)=0.0
     end if
     a(i,i)=a(i,i)+1.0_sp
  end do
  do k=n,1,-1
     do its=1,30
        do l=k,1,-1
           nm=l-1
           if ((abs(rv1(l))+anorm) == anorm) exit
           if ((abs(w(nm))+anorm) == anorm) then
              c=0.0
              s=1.0
              do i=l,k
                 f=s*rv1(i)
                 rv1(i)=c*rv1(i)
                 if ((abs(f)+anorm) == anorm) exit
                 g=w(i)
                 h=pythag(f,g)
                 w(i)=h
                 h=1.0_sp/h
                 c= (g*h)
                 s=-(f*h)
                 tempm(1:m)=a(1:m,nm)
                 a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
                 a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
              end do
              exit
           end if
        end do
        z=w(k)
        if (l == k) then
           if (z < 0.0) then
              w(k)=-z
              v(1:n,k)=-v(1:n,k)
           end if
           exit
        end if
        if (its == 30) call nrerror('svdcmp_sp: no convergence in svdcmp')
        x=w(l)
        nm=k-1
        y=w(nm)
        g=rv1(nm)
        h=rv1(k)
        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0_sp*h*y)
        g=pythag(f,1.0_sp)
        f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
        c=1.0
        s=1.0
        do j=l,nm
           i=j+1
           g=rv1(i)
           y=w(i)
           h=s*g
           g=c*g
           z=pythag(f,h)
           rv1(j)=z
           c=f/z
           s=h/z
           f= (x*c)+(g*s)
           g=-(x*s)+(g*c)
           h=y*s
           y=y*c
           tempn(1:n)=v(1:n,j)
           v(1:n,j)=v(1:n,j)*c+v(1:n,i)*s
           v(1:n,i)=-tempn(1:n)*s+v(1:n,i)*c
           z=pythag(f,h)
           w(j)=z
           if (z /= 0.0) then
              z=1.0_sp/z
              c=f*z
              s=h*z
           end if
           f= (c*g)+(s*y)
           x=-(s*g)+(c*y)
           tempm(1:m)=a(1:m,j)
           a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
           a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
        end do
        rv1(l)=0.0
        rv1(k)=f
        w(k)=x
     end do
  end do
END SUBROUTINE svdcmp_sp

SUBROUTINE svdcmp_dp(a,w,v)
  
  USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,outerprod
  USE nr, ONLY : pythag
  IMPLICIT NONE

  REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: a
  REAL(DP), DIMENSION(:), INTENT(OUT) :: w
  REAL(DP), DIMENSION(:,:), INTENT(OUT) :: v
  INTEGER(I4B) :: i,its,j,k,l,m,n,nm
  REAL(DP) :: anorm,c,f,g,h,s,scale,x,y,z
  REAL(DP), DIMENSION(size(a,1)) :: tempm
  REAL(DP), DIMENSION(size(a,2)) :: rv1,tempn

  m=size(a,1)
  n=assert_eq(size(a,2),size(v,1),size(v,2),size(w),'svdcmp_dp')
  g=0.0
  scale=0.0
  do i=1,n
     l=i+1
     rv1(i)=scale*g
     g=0.0
     scale=0.0
     if (i <= m) then
        scale=sum(abs(a(i:m,i)))
        if (scale /= 0.0) then
           a(i:m,i)=a(i:m,i)/scale
           s=dot_product(a(i:m,i),a(i:m,i))
           f=a(i,i)
           g=-sign(sqrt(s),f)
           h=f*g-s
           a(i,i)=f-g
           tempn(l:n)=matmul(a(i:m,i),a(i:m,l:n))/h
           a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
           a(i:m,i)=scale*a(i:m,i)
        end if
     end if
     w(i)=scale*g
     g=0.0
     scale=0.0
     if ((i <= m) .and. (i /= n)) then
        scale=sum(abs(a(i,l:n)))
        if (scale /= 0.0) then
           a(i,l:n)=a(i,l:n)/scale
           s=dot_product(a(i,l:n),a(i,l:n))
           f=a(i,l)
           g=-sign(sqrt(s),f)
           h=f*g-s
           a(i,l)=f-g
           rv1(l:n)=a(i,l:n)/h
           tempm(l:m)=matmul(a(l:m,l:n),a(i,l:n))
           a(l:m,l:n)=a(l:m,l:n)+outerprod(tempm(l:m),rv1(l:n))
           a(i,l:n)=scale*a(i,l:n)
        end if
     end if
  end do
  anorm=maxval(abs(w)+abs(rv1))
  do i=n,1,-1
     if (i < n) then
        if (g /= 0.0) then
           v(l:n,i)=(a(i,l:n)/a(i,l))/g
           tempn(l:n)=matmul(a(i,l:n),v(l:n,l:n))
           v(l:n,l:n)=v(l:n,l:n)+outerprod(v(l:n,i),tempn(l:n))
        end if
        v(i,l:n)=0.0
        v(l:n,i)=0.0
     end if
     v(i,i)=1.0
     g=rv1(i)
     l=i
  end do
  do i=min(m,n),1,-1
     l=i+1
     g=w(i)
     a(i,l:n)=0.0
     if (g /= 0.0) then
        g=1.0_dp/g
        tempn(l:n)=(matmul(a(l:m,i),a(l:m,l:n))/a(i,i))*g
        a(i:m,l:n)=a(i:m,l:n)+outerprod(a(i:m,i),tempn(l:n))
        a(i:m,i)=a(i:m,i)*g
     else
        a(i:m,i)=0.0
     end if
     a(i,i)=a(i,i)+1.0_dp
  end do
  do k=n,1,-1
     do its=1,30
        do l=k,1,-1
           nm=l-1
           if ((abs(rv1(l))+anorm) == anorm) exit
           if ((abs(w(nm))+anorm) == anorm) then
              c=0.0
              s=1.0
              do i=l,k
                 f=s*rv1(i)
                 rv1(i)=c*rv1(i)
                 if ((abs(f)+anorm) == anorm) exit
                 g=w(i)
                 h=pythag(f,g)
                 w(i)=h
                 h=1.0_dp/h
                 c= (g*h)
                 s=-(f*h)
                 tempm(1:m)=a(1:m,nm)
                 a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
                 a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
              end do
              exit
           end if
        end do
        z=w(k)
        if (l == k) then
           if (z < 0.0) then
              w(k)=-z
              v(1:n,k)=-v(1:n,k)
           end if
           exit
        end if
        if (its == 30) call nrerror('svdcmp_dp: no convergence in svdcmp')
        x=w(l)
        nm=k-1
        y=w(nm)
        g=rv1(nm)
        h=rv1(k)
        f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0_dp*h*y)
        g=pythag(f,1.0_dp)
        f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
        c=1.0
        s=1.0
        do j=l,nm
           i=j+1
           g=rv1(i)
           y=w(i)
           h=s*g
           g=c*g
           z=pythag(f,h)
           rv1(j)=z
           c=f/z
           s=h/z
           f= (x*c)+(g*s)
           g=-(x*s)+(g*c)
           h=y*s
           y=y*c
           tempn(1:n)=v(1:n,j)
           v(1:n,j)=v(1:n,j)*c+v(1:n,i)*s
           v(1:n,i)=-tempn(1:n)*s+v(1:n,i)*c
           z=pythag(f,h)
           w(j)=z
           if (z /= 0.0) then
              z=1.0_dp/z
              c=f*z
              s=h*z
           end if
           f= (c*g)+(s*y)
           x=-(s*g)+(c*y)
           tempm(1:m)=a(1:m,j)
           a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
           a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
        end do
        rv1(l)=0.0
        rv1(k)=f
        w(k)=x
     end do
  end do

END SUBROUTINE svdcmp_dp
